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Relativistic magnetohydrodynamics in dynamical spacetimes: 
A new AMR implementation 

Zachariah B. Etienne, Yuk Tung Liu, and Stuart L. Shapirco 
Department of Physics, University of Illinois at Urbana- Champaign, Urbana, IL 61801 

We have written and tested a new general relativistic magnetohydrodynamics (GRMHD) code, 
capable of evolving MHD fluids in dynamical spacetimes with adaptive-mesh refinement (AMR). 
Our code solves the Einstein-Maxwell-MHD system of coupled equations in full 3+1 dimensions, 
evolving the metric via the Baumgarte-Shapiro Shibata-Nakamura (BSSN) formalism and the MHD 
and magnetic induction equations via a conservative, high-resolution shock-capturing scheme. The 
induction equations are recast as an evolution equation for the magnetic vector potential, which 
exists on a grid that is staggered with respect to the hydrodynamic and metric variables. The 
divergenceless constraint V ■ B = is enforced by the curl of the vector potential. Our MHD scheme 
is fully compatible with AMR, so that fluids at AMR refinement boundaries maintain V • S = 0. In 
simulations with uniform grid spacing, our MHD scheme is numerically equivalent to a commonly 
used, staggered-mesh constrained-transport scheme. We present code validation test results, both 
in Minkowski and curved spacetimes. They include magnetized shocks, nonlinear Alfven waves, 
cylindrical explosions, cylindrical rotating disks, magnetized Bondi tests, and the collapse of a 
magnetized rotating star. Some of the more stringent tests involve black holes. We find good 
agreement between analytic and numerical solutions in these tests, and achieve convergence at the 
expected order. 
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I. INTRODUCTION 

Magnetized fluids in dynamical, strongly curved space- 
times play a central role in many systems of current inter- 
est in relativistic astrophysics. Such fluids may generate 
gamma-rays in gamma-ray bursts (GRBs), destroy difi^er- 
ential rotation in nascent neutron stars arising from stel- 
lar core collapse or binary neutron star merger, form jets 
and influence disk dynamics around black holes, affect 
magnetorotational collapse of massive stars, etc. Many of 
these systems are promising sources of gravitational radi- 
ation for detection by laser interferometers such as LIGO, 
VIRGO, TAMA, GEO and LISA. Some also emit elec- 
tromagnetic radiation, such as gamma-ray bursts, mag- 
netized disks around black holes in active galactic nu- 
clei (AGNs) and quasars, and binary supermassive black 
holes coalescing in ambient magnetized plasma. Accu- 
rate, self-consistent modeling of these systems requires 
a computational scheme capable of simultaneously ac- 
counting for magnetic flelds, relativistic magnetohydro- 
dynamics (MHD) and relativistic gravitation. 

Over the past several years, we have developed a ro- 
bust numerical scheme in 3+1 dimensions that evolves 
the Einstein equations of general relativity for the grav- 
itational held (metric), coupled to the equations of rel- 
ativistic MHD for the matter and Maxwell's equations 
for a magnetic field [l|. Our approach is based on 
the BSSN (Baumgarte-Shapiro-Shibata-Nakamura) for- 
malism to evolve the metric [2j S); ^ high- resolution. 
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shock-capturing (HRSC) scheme to handle the fluids, 
and a constrained-transport scheme to treat magnetic 
induction [J|. This GRMHD code has been subjected 
to a rigorous suite of numerical tests to check and cal- 
ibrate its validity jl|. The code has been applied to 
explore a number of important dynamical scenarios in 
relativistic astrophysics, including the collapse of magne- 
tized, differentially-rotating hypermassive neutron stars 
to black holes _5-7], the collapse of rotating stellar cores 
to neutron stars @, the collapse of rotating, supermassive 
stars and massive Population HI stars to black holes [0] , 
magnetized binary neutron star merger [lO[ , binary black 
hole- neutron stars [111, IIJIj cind the merger of binary 
black holes in gaseous environments [l3|. The purpose 
of this paper is to present a generalization of our current 
GRMHD scheme that is compatible with adaptive mesh 
refinement (AMR). 

Many problems in relativistic astrophysics require nu- 
merical simulations covering a large range of lengthscales. 
For example, to follow the final merger of a compact 
binary system with a total mass M, a lengthscale of 
^ M/30 needs to be resolved to treat the strong-field, 
near-zone regions reliably. On the other hand, accurate 
gravitational wave calculations at lengthscale ^ M must 
be performed far in the weak-field wave-zone at radius 
r > lOOM. AMR allows for sufficient resolution to be 
supplied to areas of the computational domain as needed, 
thus enabling us to resolve strong- and weak-field do- 
mains efficiently. 

One of the most subtle issues in evolving the MHD 
equations is the preservation of the divergenceless con- 
straint (V ■ B — 0) during the evolution. When evolving 
the induction equations, numerical truncation error leads 
to violations of the divergenceless constraint, resulting 



in unphysical plasma transport orthogonal to the mag- 
netic field, as well as violations of energy and momen- 
tum conservation (see e.g., 0, [iJ, [la|)- In simulations 
using a uniformly spaced grid, "constrained-transport" 
schemes (see e.g., [J, |l6|) are commonly used to main- 
tain the divergenceless constraint. In these schemes, 
special finite-differencing representations of the induc- 
tion equations are implemented to preserve a particular 
numerical representation of the divergence of the mag- 
netic field to roundoff error. In simulations using AMR 
grids, both constrained-transport schemes and the hy- 
perbolic divergence-cleaning scheme [13, Il3 have been 
used. In the hyperbolic divergence-cleaning scheme, a 
generalized Lagrange multiplier (GLM) scalar is coupled 
to the system of MHD and induction equations. No spe- 
cial finite-differencing treatment is needed in solving the 
GLM-MHD system of equations. When they appear, di- 
vergence errors of the magnetic field are both propagated 
and damped away in the GLM-MHD scheme. 

In the development of our AMR GRMHD code, we 
first tried the hyperbolic divergence-cleaning scheme, due 
to its straightforward implementation. We found that 
this scheme works well in the absence of black holes. 
One of the most commonly-adopted methods for evolving 
black holes is the moving puncture technique [19 . 20] , in 
which the physical singularity in the black hole interior 
is avoided by the use of the puncture gauge conditions. 
However, a coordinate singularity is present in the com- 
putational domain around which accurate numerical evo- 
lution is difficult to achieve. It has been demonstrated 
that the BSSN scheme, coupled with the puncture gauge 
conditions, guarantee that any inaccurate data in the 
black hole interior will not propagate out of the hori- 
zon [2l| - |23| . We find that this property is preserved in 
the presence of hydrodynamic matter. However, it is 
no longer the case in the GLM-MHD scheme. In fact, 
we find that even in the Cowling approximation in which 
the metric is fixed, inaccurate data in the black hole inte- 
rior can propagate out of the horizon in the GLM-MHD 
systems of equations. This problem may be overcome 
via black hole excision and applying appropriate ingoing 
boundary conditions at the excision boundary. (See [24 1 
for a discussion of constraint preserving boundary condi- 
tions for Newtonian MHD.) 

In developing an algorithm for maintaining V • i? = 
that is compatible with the moving puncture technique, 
we focused on constrained-transport schemes. That was 
the approach adopted in our earlier unigrid implemen- 
tation [ij. A uniform-resolution, constrained-transport 
scheme may be used on each individual AMR refinement 
level. However, maintaining the divergenceless constraint 
at refinement level boundaries requires that special inter- 
polations be performed during prolongation/restriction. 
Such pro lon g ation/restriction operators have been de- 
vised [25J, |26[, but must be fine-tuned to the particu- 
lar AMR implementation. In this paper, we propose 
an alternative, AMR-compatible constrained-transport 
scheme. Our scheme is based on the constrained- 



transport scheme described in [27[. In this scheme, 
the magnetic induction equation is recast as an evolu- 
tion equation for the magnetic vector potential. The 
divergence-free magnetic field is computed via the curl of 
the vector potential. The evolution of the vector poten- 
tial is carried out in the same HRSC framework as other 
hydrodynamic variables. This scheme is numerically 
equivalent to the commonly used constrained-transport 
schemes based on a staggered mesh algorithm [lg|. This 
scheme is readily generalized to an AMR grid. Unlike the 
magnetic field, the vector potential is not constrained, 
and so any interpolation scheme can be used during pro- 
longation and restriction, thus enabling its use with any 
AMR algorithm. 

We have performed several tests on our new AMR 
constrained-transport scheme. We find that it works well 
even in black-hole spacetimes. Inaccurate data generated 
in the black hole interior stay inside the horizon. Hence 
our scheme is compatible with the moving puncture tech- 
nique. 

The structure of this paper is as follows. In Sec. [Hi 
we describe our formalism, focusing on the derivation of 
the evolution equation for the magnetic vector poten- 
tial. Then we describe our numerical scheme to evolve 
the coupled Einstein-Maxwell-MHD equations fSec.lIIip. 
Next we present several stringent code tests, including 
one- and two-dimensional shocks, magnetized Bondi ac- 
cretion and the collapse of a magnetized rotating star 
fSec. lIV| ). Finally, we summarize our work in SecFVland 
discuss applications of our new code to study various in- 
teresting problems in relativistic astrophysics. 



II. FORMALISM 

The formulation and numerical approach adopted in 
this paper are basically the same as those already re- 
ported in our previous work [l], [Tl], [I^], to which the 
reader may refer for details. Here we introduce our nota- 
tion, summarize our method, and focus on the derivation 
of the evolution equation for the magnetic vector poten- 
tial in the ideal MHD limit, which is the basis of our new 
AMR constrained-transport scheme. Geometrized units 
{G ~ c — I) are adopted throughout. Greek indices de- 
note all four spacetime dimensions (0, 1, 2, and 3), and 
Latin indices imply spatial parts only (1, 2, and 3). 



A. Metric evolution and gauge conditions 

We use the standard 3-1-1 formulation of general rela- 
tivity and decompose the metric into the following form: 



ds' 



-a^df + j,j{dx' + p'dt)(dx^ -f p'dt) 



(1) 



The fundamental variables for the metric evolution are 
the spatial three-metric ^ij and extrinsic curvature Kij. 
We adopt the BSSN formalism [3, Q in which the evolu- 
tion variables are the conformal exponent (f> = ln(7)/12. 



the conformal 3-metric 7^ = e '^'^jij, three auxihary 
functions F' = — 7*"' ,j, the trace of the extrinsic curva- 
ture K = jijK'^^, and the trace- free part of the con- 
formal extrinsic curvature Aij = e~'^'^{Kij — "fijK/?)). 
Here, 7 = det(7ij) is the determinant of the spatial 
metric. The full spacetime metric g^n, is related to the 
three-metric 7^1/ by 7^1/ = g^y +n^nv, where the future- 
directed, timelike unit vector n'^ normal to the time slice 
can be written in terms of the lapse a and shift /3* as 
n^^ — q;^^(1, — /3*). The evolution equations of these 
BSSN variables are given by Eqs. (9)-(13) in [il|. It 
has been suggested that evolving x = e"^"^ or VF = e"^*^ 
instead of 4> gives more accurate results in binary black 
hole simulations (see e.g. [28-30]). Our code is capa- 
ble of evolving these variables. Kreiss-Oliger dissipation 
is sometimes added in the BSSN evolution equations to 
reduce high-frequency numerical noise associated with 
AMR refinement interfaces [3l|. It is also found that 
Kreiss-Oliger dissipation is sometimes useful in hydrody- 
namic simulations involving a black hole in a dynamical 
spacetime |12l.l32|. 

We adopt standard puncture gauge conditions to 
evolve the lapse and shift: an advective "l-f-log" slicing 
condition for the lapse and a "Gamma-freezing" condi- 
tion for the shift [33|. The evolution e qua tions for these 
quantities are given by Eqs. (2)-(4) in [l3 |. 



B. Evolution of electromagnetic fields 



by 



The electromagnetic stress-energy tensor TJ^ is given 



where the specific enthalpy h is related to the specific 
internal energy ehy h = \ -\- e + P/ Po- The electric and 
magnetic fields measured by an observer comoving with 
the fluid are [cf. Eq. (g])] 



(u) " 



B 



(u) 



p*I^tJ, 



(7) 



For many applications of interest in relativistic astro- 
physics, one can assume perfect conductivity. In this 
ideal MHD limit. Ohm's law yields the MHD condition: 



.F^" = , 



(8) 



which is equivalent to the statement that the electric field 
observed in the fluid's rest frame vanishes {Et^. = 0). In 
this limit, the total stress-energy tensor is given by 

&2 



T^"" = {p^h + b^)ui'u'' + P + 



2 



where 6'^ = B'^./\/At: and h — u w^^ 
related to B^^ by (see [1] for a derivation) 



The vector h^^ is 



6^' = -- 



Pf'^B" 



TlyU'^^/A'IT 



(10) 



where P^^ = g^^ + u^u^ is a projection tensor. 

The evolution equation for the magnetic field in a per- 
fectly conducting MHD fluid can be obtained in conser- 
vative form by taking the dual of Maxwell's equation 



F, 



[/ji/.A] 



0. One finds 



y^p*f^i^ 



-.dAv^F*'^n^o, 



(11) 



Te'Z = ^ [f'^'F'^^ ~ \g^''F^pF'^P ) . (2) 

We decompose the Faraday tensor F^'' as 

F^" = n^F"^ - n^'E^' + n^e'^'^^Bs , (3) 

so that E^ and B^ are the electric and magnetic fields 
measured by an observer normal to the spatial slice n^. 
Both fields are purely spatial [E^^n^ = B^Uf^ = 0), and 
one can easily show that 

F^ = F^-n, , B^ = ^e^--^n,F^, = n.F*'''' , (4) 



where 



p*t2iy _ Zft^^i^^E 

2 



(5) 



is the dual of F^"' . 

Along with the electromagnetic field, we also assume 
the presence of a perfect fluid with rest-mass density poi 
pressure P, and 4- velocity u^, so that the total stress- 
energy tensor is 



where ^/— 5 = ot^Py. The time component of Eq. (jlip 
gives the no-monopole constraint 



where 



djP = , 



P' - ^B' 



(12) 



(13) 



The spatial components of Eq. ([Tl]) give the magnetic 
induction equation, which can be written as 



where v^ — u^ /v!^ . 

The induction equation can be recast as 

dtB'^l'^'^himdMBn > 



(14) 



(15) 



where both e'-''^ and eijk denote the permutation symbol, 
i.e. they are equal to 1 if ijk arc in even permutation 
of (1,2,3), —1 if in odd permutation, and if any two of 
the indices are equal. The divergenceless constraint p^ 
implies that P* can be derived from a vector potential 
A,:. 



T^'^ = pohu^'u'' + Pgt''' + TZ , 



(6) 



B' 



?^"9,^fe 



(16) 



It follows from Eqs. ^3^ and ([16]) that 



B' 



f'^'^d.Ak 



(17) 



n^^m^ is the three-dimensional 



where e^^^ = l^^^ j ^ 
Levi-Civita tensor associated with 7^. Equation ()17p can 
be derived in a more general framework, as shown in [3J] 
The induction equation ([15]) will be satisfied automat 
ically if Ai satisfies the evolution equation 



dtAi = hjkV^B^ 



(18) 



It is clear that the evolution equations for Ai are not 
unique, since there are gauge degrees of freedom in the 
electromagnetic 4-vector potential. The general evolu- 
tion equation for Ai in the ideal MHD limit is obtained 
by combining Eqs. (33) and (46) in [34]: 



dtA^ = h.kv'B'^ - a.(a$ - ^'A,) 



(19) 



where <I> is the electromagnetic scalar potential. Hence 
the evolution equation (|18l) is equivalent to choosing the 
electromagnetic gauge condition 






(20) 



where C is a constant. In Minkowski spacetime, in which 
a = 1 and /3* = 0, the gauge condition reduces to <I> = C 

The scalar potential is only needed if one wishes to 
compute the electric field E^ . However, in the ideal MHD 
limit, the condition Ufj^F'^" = relates E^ to B' and v'^: 
aEi = —£ijk{v^ + (i-')B^ . Therefore, it is not necessary 
to keep track of the scalar potential $ in the ideal MHD 
limit. 

In the nonrelativistic limit, Eq. ([TS]) reduces to 



dtB^V x{vx B) 
and Eqs. ([T7]) and ([18]) reduce to 



B = V xA 



dtA ~ V X B 



(21) 



(22) 



The ideal MHD condition becomes E — —v x B. 

In our new AMR constrained-transport scheme, the in- 
duction equation is evolved via Eq. (|18p . The divergence- 
free magnetic field is then computed using Eq. ([T7)) . The 
numerical implementation will be described in Sec. IIIII 



TABLE L Storage location on grid of the magnetic field B^ 
and vector potential Ai 

Variable storage location 



B^B^ (i + ii,fe) 

By, By (i,j + i,fc) 

B\B^ (i,j,fc + i) 



Ay 

A, 



(i+j.J + jfc) 



The evolution equations are derived from the rest-mass 
conservation law ^^(pqu^) = and conservation of 
energy-momentum V^T'*" = 0. These result in the con- 
tinuity, momentum and energy equations (Ij 



dtp* + dj{p*iP) = , 



(27) 



dtS, + d,ia^Th) = -aV7r"^5a/3,. , (28) 

dtf + d,{a^y^T"'-p*v') = s, (29) 

where the source term in the energy equation is given by 

= a^l [(r"°/3*/3^' + 2r"'^^' + T'^)K,j 

-(r™/3* + r°')a,a] . (so) 

To complete the system of equations, the fiuid equation of 
state (EOS) is specified. Our code currently implements 
a hybrid EOS of the form [35,] 

P{Po, e) = -Pcoid(po) + (rth - l)po[e - ecoid(po)] , (31) 

where Pcoid and ecoid denote the cold component of P 
and e respectively, and Tth is a constant parameter which 
determines the conversion efficiency of kinetic to thermal 
energy at shocks. The function ecoid(po) is related to 
Pcoid(po) by the first law of thermodynamics. 



, X 1 -Pcold(Po) , 

ecoid(po) = / 2 ^Po 

Po 



(32) 



In the code tests presented in this paper, we adopt the 
F-law EOS P = (r — l)poe. This corresponds to setting 
Pcoid = i^Pq (with constant k) and Fth = T. 



C. Evolution of the hydrodynamics fields 

The stress-energy tensor for a magnetized plasma in 
the ideal MHD limit is 

Tf"" = {poh + b^juf^u" +(P+J) 9'"" - b^b" . (23) 



Our evolution variables are 



P* = -y/lPonf.u'', 



(24) 
(25) 
(26) 



III. NUMERICAL IMPLEMENTATION 

We adopt Cartesian coordinates in our 3-1-1 simula- 
tions. Equatorial symmetry (i.e. symmetry with respect 
to the reflection z -^ —z) is imposed when appropri- 
ate to save computational time. All the BSSN and hy- 
drodynamical variables are stored at grid points (i, j, k). 
Magnetic field B^ and vector potential Ai are stored at 
staggered grid points as summarized in Table HI 

The BSSN equations are evolved using a finite- 
differencing scheme. Our code currently supports sec- 
ond, fourth, and sixth order spatial finite-differencing. 



In a spacetime containing black holes, we typically use 
a fourth or sixth order finite-differencing scheme. Our 
code is embedded in the Cactus parallelization frame- 
work [3a |. with time-stepping managed by the MoL 
(Method of Lines) thorn, which supports various explicit 
time-stepping algorithms. Typically, we use the fourth- 
order Runge-Kutta method in time when evolving space- 
times containing black holes. 

We use the Carpet [33| infrastructure to implement 
moving-box adaptive mesh refinement. In all AMR sim- 
ulations presented here, second-order temporal prolon- 
gation is employed, coupled with fifth-order spatial pro- 
longation for evolution variables stored on the unstag- 
gered grid. The memory allocation for the staggered 
variables are the same as the unstaggered ones. The 
staggering is incorporated in our code in the evolu- 
tion steps. Different spatial prolongation and restric- 
tion schemes have to be applied on the staggered evo- 
lution variables Ai to account for the different relative 
positions of these variables on adjacent refinement lev- 
els. We currently use a third-order Lagrangian scheme 
for interpolating these variables, but it can be eas- 
ily generalized to other higher-order schemes, as well 
as more sophisticated schemes such as the essentially 
non-oscillatory fENC) |38| and weighted essentially non- 
oscillatory (WENO) ^^ |40] schemes. We plan to inves- 
tigate these alternative schemes in the future. 



conservative variables U from the point-value primitive 
variables P using a finite-volume algorithm. Next the 
updated point-value U is computed from the updated 
volume average quantity U to the desired order of ac- 
curacy. The updated point-value P is then computed 
from the updated point-value U and metric quantities 
through primitives inversion. In this paper, we only con- 
sider second-order schemes for simplicity. Higher-order 
schemes are planned for the future, and important ex- 
tra steps necessary to go beyond second-order will be 
reviewed in this section. 

Equation (|33p can be written in a finite-volume form 
by integrating it over a cell volume. We obtain 



dtU 



fJ i,j,k 



(A.(F)),,,-fe , (A,(F)),,,-fe 



Aa; 

(A. (F) ),,,., 






where 



(a,(f)),,;,^(f; 



i+ij,fc 



i-:jj",fc 



(34) 



(35) 



and similarly for operators Aj^ and A^. We note that 
only a subset of U, i.e. p*, f and Si, is evolved using 
Eq. (pi)) . The evolution of B^ will be described in the 
next subsection. The bracket () denotes a surface aver- 
age. For example. 



1 



A. MHD evolution 

The technique for evolving the BSSN equations is de- 
scribed in our earlier papers [ll|, [Ij, |41|, so we focus 
here on our MHD evolution technique, which is based 
on an HRSC scheme. The goal of this part of the nu- 
merical evolution is to determine the fundamental MHD 
variables P = {pQ,P,v'^,B'^), called the "primitive" vari- 
ables, at future times, given initial values of P. The 
evolution equations ([H]), (P7)) -(pn | are written in con- 
servative form: 



dtU + V ■F = S , 



(33) 



where U{P) = (p*,f, S'i, J3*) are the "conserved" vari- 
ables, and the flux F{P) and source S{P) do not contain 
derivatives of the primitive variables, although they are 
explicit functions of the metric and its derivatives. 

Equation (j33[) may be evolved using a flnite- volume or 
finite-difference scheme. A finite-volume scheme evolves 
the volume-averaged variables, whereas a finite-difference 
scheme evolves the point-valued variables. Our adopted 
constrained-transport scheme is based on a finite- volume 
algorithm. In a second-order scheme, there is no distinc- 
tion between these two types of methods since the volume 
average and the gridpoint value are the same to second 
order. Since the metric is evolved using a finite-difference 
scheme, care must be taken to evolve the MHD and 
induction equations using a higher-order finite-volume 
scheme. One solution is to evolve the volume averaged 



^+^.J,k 



AyAz 



dy 



dzF 



xi , y, z 



(36) 



where xf = Xi ± Aa;/2, yf = yj ± Ay/2 and z^ = Zk ± 
Az/2. The fluxes {F)^j_^_ij. and (F)jj j,_|_i are defined in 
the same way except that the surfaces to be averaged are 
in the x-z plane and x-y plane, respectively. The surface 
averaged fiux (F) and the point- value flux F are the same 
to second-order accuracy. To implement a higher-order 
scheme, one needs to compute not only the point- value 
F at the zone interface to the desired order, but also (F) 
from the point- value F to the desired order of accuracy. 
The computation of the fluxes is basically the same 
as described in [1]. It involves the reconstruction step 
and the Riemann solver step. In the reconstruction 
step, primitive variables at the zone interface are re- 
constructed. A slope-limited interpolation scheme from 
the zone center gives Pr and Pl, the primitive vari- 
ables at the right and left side of each zone interface, 
respectively. We usually employ the piecewise parabolic 
method (PPM) [il] or the monotonized central (MC) \^ 
reconstruction scheme, but in some problems involving 
strong discontinuities a more diffusive scheme such as 
the minmod reconstruction scheme must be used (see 
Sec. lIVA"2|) . Since B^ is staggered (as shown in Table IJ), 
each B* at one of the zone interfaces need not be com- 
puted. From Pfl and Pl , we compute the fluxes Fji and 
Fl, the "conservative" variables Ur and Ul, as well as 
two pairs of characteristic velocities c^ and c±. at each 
zone interface (see Sec. IIIB of [1| for details). 



The next step is the Riemann solver step. We em- 
ploy the HLL (Harten, Lax, and van Leer) approximate 
Riemann solver [4J] in which the HLL flux is given by 



pHLL _ 



c-Fn 



Fl - C+C {ur " Ul) 



(37) 



where c^ = max(0, ±c^, ±c^). Our code also has the op- 
tion of using the single-speed, local Lax-Friedrichs (LLF) , 
or central-upwind, flux. 



F 



LLF 



1 



[Fr + Fl-c(ur-ul)] 



(38) 



where c 

The accuracy of the resulting flux depends on the re- 
construction scheme and Riemann solver. In a smooth 
flow, MC reconstruction results in a second-order ac- 
curate point-value flux F, whereas PPM is third-order. 
However, these two schemes reduce to first-order in a 
discontinuous flow (e.g. shocks) or at local extrema of 
P. As mentioned above, even in a smooth flow where 
PPM gives third-order accurate point-value F, (F) has 
to be computed from the point-value F to third order to 
achieve an overall third-order accuracy. 



B. Constrained transport scheme 

In this subsection, the standard constrained-transport 
scheme based on the staggered algorithm [16| is reviewed 
briefly. Next we introduce the vector potential method 
described in [231 ■ These two approaches give numerically 
identical results for schemes in which the time integration 
and spatial derivatives commute. 

The evolution variables for the magnetic field in the 
standard constrained-transport scheme is the surface av- 
eraged field {B^) defined in the same way as the surface 
averaged fluxes: 






1 r^ 



„+ „,+ 



1 



■2''^ AxAz 



dy dzB-{xt,y,z)i39) 

dx dzBy{x,y+,z)(4Q) 

dx [ ' dyB^x,y,z+%Al) 



X, Jy 



Integrating the magnetic constraint equation djB^ = 
over a cell volume gives the finite- volume equation for 
the constraint 



(A,(S-)),,,-fe , iAy{By)),,j,k , (A,(^-)),,,-fe 



Ax 



Ay 



Az 



(42) 



To derive the finite- volume equation for the magnetic in- 
duction equation, we first rewrite Eq. (|15p as 



where 



%B^ = 


-dy£' + d^£y , 


(43) 


%By = 


-d,£' + d,£' , 


(44) 


%B' = 


-d^sy + dyE^ , 


(45) 


f" = 


-vyJB^ + v^By , 


(46) 


£y = 


-v'-B^ +v''B' , 


(47) 


£' = 


-v'^By + vyB"" . 


(48) 



We next define the line averaged £^ as 



ex 



£y 






^/ f-(x,y+,z+)da:, (49) 



1 

1 

Ai 



f^(x+,2/,z+)dy, (50) 
£-(a;+,y+,z)dz. (51) 



Note that f * is staggered in the same way as Ai (see Ta- 
ble H]). The finite- volume equations for the magnetic in- 
duction are obtained by integrating Eq. (P5| over the cell 
surface normal to the x-direction, integrating Eq. ((44]) 
over the cell surface normal to the y-direction, and in- 
tegrating Eq. (HS)) over the cell surface normal to the 
z-direction: 






(A,g^),+ .,^.fe (A,g-),^l^^.fc 

Az Ay 

(A.g-),^^-+.,fc (A,g-),,,^.+ i^, 

Ax Az 

(A,5-),,^.,+ i {Ady\^,^^. 



Ay 



,(52) 
,(53) 
.(54) 



It is straightforward to verify that Eqs. ((52|) - ([54| im- 
ply that the time derivative of the left hand side of 
Eq. P^ vanishes. Hence a finite- volume scheme that 
evolves Eqs. ([5^ -([M | preserves the constraint P^ to 
roundoff error, provided that the initial data (i?*) satisfy 
the constraint. 

To evolve Eqs. dS^ -lfM ]) . £ has to be computed at 
the zone edge. The computation is similar to that of 
the flux F described in the previous subsection. Since 
B^ is staggered (as specified in Table HI, computation 
of each £* at the zone edge requires reconstruction of 
B^ along one direction. However, since v^ is stored at 
the zone center, two independent one-dimensional recon- 
structions are necessary, as pointed out in ^] . The HLL 
and Lax-Friedrichs formulas for £^ at the zone edge are 
given by (27 1 



r) 



HLL 



'^x ^y ^LL 



^x Cy £lR 



^x Cy £rl 



'^x S £rr 



(Ct + Cx )iCy + Cy ) 



Cx ~r Cx 



{Bl 



and 



in 



z\LLF 



. C, 



Bl) 



■■LR 



c c 



■{B^R 



Bl) 



(55) 



The evolution equation is derived from Eq. 
given by 



181) and is 



-TVL 



^rr) 



^{Bl-BD-^-^iBf 



Bl) , (56) 



which are the generahzations of Eqs. (|37)) and (|38)). In 
the above formulas, f£^ denotes the reconstructed left 
state in the x-direction and right state in the y-direction. 
Other symbols involving £^ are interpreted in the similar 
fashion. B^ and Bj^ denote the reconstructed right and 
left state of B^ in the x-direction; Bf^ and B£ denote the 
reconstructed right and left state in the y-direction. The 
c^ and Cy should be computed by taking the maximum 
characteristic speed among the four reconstructed states. 
However, we set them equal to the maximum over the two 
neighboring interface values for simplicity, as suggested 
in |27| . In the LLF formula, Cx and Cy are set to the 
maximum of c^ and c^, respectively. The formula for 
^^K-jHLL jg obtained from Eq. (ISS)) by permuting the in- 
dices z -^ X, X ^ y and y ~> z, whereas the formula 
for {£y)^^^ is obtained from Eq. (|55p by permuting the 
indices z — )■ y, a; — >■ z and y ^ x. The same rule applies 
for (f ^)LLF and {£y)^^^ . The reconstructed point-value 
£^ at the zone edge is the same as the line averaged value 
f to second-order. If one wishes to go beyond second- 
order, £' has to be computed from £^ to the desired order 
of accuracy. 

We now describe the vector potential method proposed 
in [23| , which has been adopted for our AMR constrained- 
transport scheme. We first define the line averaged vector 
potential Ai exactly the same way as £'': 



(^.).. 



1 



j+i fc+i 



Axix,y+,z^)dx , (57) 



1 /"^j 

1 f^^ 
(iz),+ i,,+ i,fe = ^ I AAx+,y+,z)dz . (59) 



It follows from Eq. dTH) that 

~ (Ayi.).+|,j-fc (A.4).+ i,j-,fc 
[B )i+^j^k = TTT. TZ ' i^^) 



(By) 

(B 



Ay Az 

''^■+5''^- ~ Az Ax 



i,J,k+i 



Ax 



Ay 



,(61) 
.(62) 



It is easy to verify that the data (_B*) generated from Ai 
from the above formulas satisfy the constraint (|^^ . In 
the vector potential method, the evolution variable is Ai. 



5t(^a;)»j+i,fc+i = 


cx 


(63) 


5t(iy)j+i,j-fc+i = 


cV 


(64) 


9t(A^)j+i,j+i,fe = 




(65) 



The value of £^ at the zone edge is computed in ex- 
actly the same way as the standard constrained-transport 
scheme, i.e., by using Eq. (1551) or Eq. ((55|) for £' and 
similar formulas for £^ and £y . Having evolved Ai, the 
updated conservative variables (B*) are computed using 
Eqs. (|5D l) - (p^ . The divergence of {B^) is therefore auto- 
matically guaranteed to be zero to roundoff error. 

It is apparent that the vector potential method and 
the standard constrained-transport scheme are closely re- 
lated. They both apply the same procedure of computing 
f at the zone edge. They both involve taking spatial 
derivatives (more precisely, the discretized curl opera- 
tor) via the differencing operators Aj;, Aj, and A^. The 
only difference between these two methods is that in the 
standard constrained-transport scheme, spatial deriva- 
tives are applied before time integration, whereas in the 
vector potential method spatial derivatives are applied 
after time integration. Since we employ the MoL algo- 
rithm in which spatial derivatives and time integration 
commute, the two methods give numerically identical re- 
sults in simulations using a uniformly-spaced grid. We 
prefer to use the vector potential method in AMR sim- 
ulations since Ai is not constrained and so does not re- 
quire special interpolation schemes during prolongation 
and restriction. 

During the MHD evolution steps, values of B^ at the 
zone center are also needed, which are currently com- 
puted by simply taking the average of B^ on the stag- 
gered grid. Taking the limit Ax' — > in Eq. (|42p . we 
see that (B*) is always continuous in the x' direction 
(e.g., even in the presence of shocks). Thus the averag- 
ing scheme generally gives a second-order accurate B^ at 
the zone center. Higher-order schemes will require more 
sophisticated interpolation algorithms. 



C. Recovery of primitive variables 

Having computed U at the new timestep, we need to 
recover P, the primitive variables on the new time level. 
This is not trivial because, although the relations U(P) 
are analytic, the inverse relations P(U) are not. For a 
F-law EOS P — (T — l)poe- the inversion can be reduced 
to an eighth-order polynomial equation [23, |43 . In the 
absence of magnetic field, the equation can be further 
reduced to a quartic equation where an analytic solution 
is available. However, for a general EOS, the inversion 
must be solved numerically. Various inversion algorithms 
are studied extensively in [45|, and it has been found 
that the most efficient inversion technique is to solve two 
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coupled nonlinear equations using the Newton-Raphson 
scheme. 

Our code supports three inversion schemes: the op- 
timal 2D scheme described in [451 , a slightly modified 
analytic quartic solver from the GNU Scientific Library 
(used for a F-law EOS in the absence of magnetic field) , 
and our older scheme that solves four coupled nonlinear 
algebraic equations. 



a Kreiss-Oliger dissipation in the black hole interior is 
found to be useful [ij, [S^j , as well as setting an upper 
and lower limit on the pressure. In some MHD simu- 
lations, we find that setting the magnetic field to zero 
deep inside the horizon can stabilize the evolution (see 
Sec. WUi . 



E. Low-density regions 



D. Black hole interior 

We use the moving puncture technique to evolve space- 
times containing black holes. The black hole spacetime 
singularity is avoided by the puncture gauge conditions, 
but a coordinate singularity (i.e. puncture) remains in 
the interior of each black hole on the computational do- 
main. One nice property of the moving puncture method 
is that, although accurate evolution near the puncture 
is not maintained, inaccurate data do not propagate 
out of the black hole horizon. This method proves to 
be robust in the evolution of binary black holes and is 
widely used in the numerical relativity community. The 
moving puncture method has also been used in simula- 
tions involving hydrodynamic matter and MHD (see e.g., 

EmEllall&Si). 

One difficulty in handling MHD in the black hole in- 
terior is the loss of accuracy near the puncture. This 
can drive the "conservative" variables U out of physical 
range, resulting in unphysical primitive variables after 
inversion (e.g. negative pressure or even complex solu- 
tions). In the absence of magnetic fields, this can be 
avoided by enforcing the constraints |llj] 

\S\^ = Y^S^Sj < f{T + 2p^),and (66) 

f > . (67) 

When the second condition is not met, we reset r to 
a small positive number. When the first condition is 
violated we rescale Si so that its new magnitude is 
l'5'P — fT{f + 2p^,), where / < 1 is a parameter. This 
technique does not apply in the presence of magnetic 
fields. We instead apply a fix, first suggested by Font 
et al [401, which consists of replacing the energy equa- 
tion (pg|) by the cold EOS, P = Pcoid(po) when solving 
the system of equations. This substitution guarantees a 
positive pressure. In rare cases, this revised system also 
fails to give a solution and we repair the zone by aver- 
aging from nearby zones (averaging is not applied to the 
magnetic field). 

When matter and magnetic fields fall into the black 
hole, the energy density near the puncture can be very 
high. This results in a large energy source term in the 
BSSN equation, which can cause the conformal related 
metric 7^ to lose positive definiteness near the puncture. 
This behavior eventually causes the code to crash. Hence, 
other techniques are sometimes used to stabilize the evo- 
lution in the black hole interior. For example, adding 



If a pure vacuum were to exist anywhere in our com- 
putational domain, the MHD approximation would not 
apply in this region, and the vacuum Maxwell equations 
would need to be solved there. In many astrophysical 
scenarios, however, a sufficiently dense, ionized plasma 
will exist outside the stars or disks, where MHD will 
remain valid in its force- free limit. As in many hy- 
drodynamic and MHD simulations, we add a tenuous 
"atmosphere" to cover the computational grid outside 
the star or disk. We maintain a density and pressure 
floor yOatm and Patm in the atmosphere. We usually set 

Patm = 10"^°Pmax(0) and Patm = -Pcold(/3atm) , where 

Pmax(O) is the maximum rest-mass density in the ini- 
tial data. Throughout the evolution, we impose limits on 
the atmospheric pressure to prevent spurious heating and 
negative values of the internal energy when the density 
Pq is smaller than a threshold pth- Specifically, we require 

■Pmin(/Oo) < P < PmnxiPo), where Pmax(po) = lOPcold(po) 

and Pmin(po) = Pcoid(/Oo)/2 when po < pth- The value of 
Pth is usually set between lOpatm and lOOpatm- Setting 
Pth too high could cause unphysical effects in a simula- 
tion, such as spurious angular momentum loss jl2|. 



F. Boundary conditions 

In simulations in which the spacetime is asymptoti- 
cally fiat, we apply Sommerfeld outgoing wave boundary 
conditions to the BSSN and gauge variables f , i.e.. 



f(r,i) 



Ar 



f (r -Ar,t- AT) 



(68) 



at the outer boundary of our numerical grid. Here AT 
is the timestep and Ar = Q;e~^'^AT. In simulations in 
which hydrodynamic matter and plasma are initially lo- 
calized, outflow boundary conditions are imposed on the 
hydrodynamic variables po, v^ and P (i.e., the variables 
are copied along the grid directions with the condition 
that the velocities be positive or zero in the outer grid 
zones). For the magnetic field, we compute Ai at the 
outer boundaries by either linear or quadratic extrapo- 
lation. The linear extrapolation is equivalent to copying 
B^ to the outer boundary, whereas the quadratic extrap- 
olation corresponds to linearly extrapolating S' to the 
outer boundary. 

Other boundary conditions are used in the code tests 
presented in this paper, which will be specified in each 
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FIG. 1: ID fast and slow shock density and velocity profiles, 
at i = tflnai (see Table |ll|. Data from numerical simulations 
with resolution Aa; = 0.01 are plotted with dashed (black) 
lines, and solid (red) lines denote the analytic solutions. 



FIG. 2: Same as Fig.[T]but for the ID switch-off and switch-on 
tests. 



IV. CODE TESTS 
A. Minkowski spacetime MHD tests 

1. One-dimensional tests 

We perform a suite of one-dimensional MHD tests in 
Minkowski spacetime, as described in [51|. The initial 
configurations of the tests are summarized in Table HH 
We only perform tests in which analytic solutions are 
available. In these ID tests, all variables are functions of 
X only. The divergenceless constraint V • B = implies 
that B^ is a constant. For these tests, the initial magnetic 
field B can be derived from the following vector potential: 



A^{x) = 



Ay{x) 



B'{x')dx' , 



A.{x) = vB- 



By{x')dx' . 



(69) 
(70) 

(71) 



We integrate the MHD equations from t = to 
i — ifinaij where ignai is specified in Table |ll] for each 
case. The gas satisfies a F-law EOS with F — 4/3, and 
is evolved on a uniform resolution grid with Aa; — 0.01. 
We are able to integrate all the cases using the PPM 
reconstruction scheme and the HLL approximate Rie- 
mann solver with a timestep Ai = O.SAx. We use 
"copy" boundary conditions (i.e. hydrodynamic variables 



are copied and the vector potential is linearly extrapo- 
lated to the boundary points) in all cases. The first 6 
tests in Table [H] start with discontinuous initial data at 
X = 0, with homogeneous profiles on either side. Fig- 
ures [TH3] show the profiles of po and u^ at time t = ignai- 
Notice that the numerical results agree very well with the 
analytic solution in all cases. The overall performance of 
the new MHD scheme in these tests is about as good as 
our old constrained-transport scheme presented in [ij. 

In the nonlinear Alfven wave test, unlike the first 6 
tests, the two states listed in Table |lT] are joined by a 
continuous function. We use the same initial data de- 
scribed in Appendix B of jlj. Figure |4] demonstrates very 
good agreement between numerical results and the ana- 
lytic solution for velocity and magnetic field profiles at 
time t — ifinai — 2. Figure [5] shows the L2 norm of the 
error in u^ , u^, B^ and B^' , varying only numerical res- 
olution. The L2 norm of a grid function Sg = g — g'^^^'^'^ 
is computed by summing over every grid point i: 



L2{Sg) 



\ 



N 



Y.[5g{x,)?^^ 



(72) 



where N ex 1/Aa; is the number of grid points. We find 
that the errors converge to zero at second order in Ax, 
as expected. 



TABLE II: Initial states for ID MHD tests; 



Test 



Left state 



Rigfit State 



ifln 



Fast Sliock u' = (25.0, 0.0, 0.0) 

ByV^= (20.0,25.02,0.0) 
P = 1.0, po = 1.0 



u' = (1.091,0.3923,0.00) 2.5 
ByV^= (20.0,49.0,0.0) 
P = 367.5, po = 25.48 



Slow Shock u' = (1.53, 0.0, 0.0) 

B7^/4^= (10.0,18.28,0.0) 
P = 10.0, po = 1.0 



m' = (0.9571, -0.6822, 0.00) 2.0 
B'/V^= (10.0,14.49,0.0) 
P = 55.36, po = 3.323 



Switch-off Fast 
Rarefaction 



u' = (-2.0,0.0,0.0) 

B7^/4^= (2.0,0.0,0.0) 

P=1.0, po = 0.1 



m' = (-0.212,-0.590,0.0) 1.0 
ByVI^= (2.0,4.71,0.0) 
P = 10.0, po = 0.562 



Switch-on Slow 
Rarefaction 



u' = (-0.765,-1.386,0.0) 
By^/I^= (1.0,1.022,0.0) 
P = 0.1, po = 1.78 X 10^=* 



u' = (0.0,0.0,0.0) 

ByV^= (1.0,0.0,0.0) 

P= 1.0, po = 0.01 



Shock Tube 1 



«' = (0.0,0.0,0.0) 

57^4^= (1.0,0.0,0.0) 

P = 1000.0, po = 1.0 



Shock Tube 2 



Nonlinear Alfven wave 



u' = (0.0,0.0,0.0) 

B7^4^= (0.0,20.0,0.0) 

P = 30.0, po = 1.0 



u' = (0.0,0.0,0.0) 

By^/I^= (3.0,3.0,0.0) 

P=1.0, po = 1.0 



2.0 



m' = (0.0,0.0,0.0) 1.0 

57^4^= (1.0,0.0,0.0) 
P = 1.0, po = 0.1 



m' = (0.0,0.0,0.0) 1.0 

57^4^= (0.0,0.0,0.0) 
P = 1.0, po = 0.1 



u* = (3.70,5.76,0.00) 2.0 

P7^/4^ = (3.0, -6.857, 0.0) 
P=1.0, po = 1.0 



'^ In all cases, the gas satisfies the F-law EOS with F = 4/3. For the first 6 tests, the 
left state refers to a; < and the right state, a; > 0. 

^ For the nonlinear Alfven wave, the left and right states are joined by a continuous 
function. See [53] or Appendix B of [l| for details. 
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2. Two-dimensional tests 

We perform the two-dimensional cylindrical blast ex- 
plosion test and rotating disk test described in [1^, [23] . 
In both tests, all variables are functions of a; and y only, 
the initial magnetic field is uniform and oriented along 
the x-direction, and the initial velocity does not have 
the z-component. Such a uniform magnetic field can be 
derived from the vector potential 



A,^Ay=0 , A,=yB^ 



(73) 



It can be shown from the MHD evolution equations 
that Ax — Ay ^ B^ = v^ ~ remains true for all time t. 
It can also been shown from the finite- volume equations 
that our MHD evolution scheme preserves this property. 
Our numerical simulations also confirm that A^ — Ay ^ 
B^ — v^ = is satisfied to roundoff error at all times. It 
follows from A^ = Ay = B^ = and B = V x A that 
B'^diAz — 0. Hence contours of constant A^ coincide 
with the magnetic field lines. The evolution of magnetic 
field thus reduces to the evolution of A^, which can be 
shown to satisfy the simple advection equation: 



dtA, -f v'd.A, 







(74) 



It follows from Eq. (fM)) that the constant Az contours 
are comoving with the fluid. We note that we do not 
evolve Eq. ((7^ directly. Instead, we evolve Az using 



the HRSC scheme described in Sec. IIIIBI Small numeri- 
cal resistivity inherent in our HRSC scheme could cause 
small violations of Eq. (|74p , especially in regions of strong 
shocks. However, the relation B'^diAz = is satisfied to 
truncation error at all times. 

Cylindrical blast explosion 

In this test, the fluid is initially at rest with uniform 
density po = 1 throughout the computational domain 
X G (-0.55,0.55), y G (-0.55,0.55). Inside a cylin- 
der of radius 0.08 is a uniform high pressure P = lO''^ 
surrounded by an ambient fluid of much lower pressure 
P = 0.1. The initial magnetic field is B^/v^ = 4.0, 
and B^ = B^ = everywhere. The fluid satisfles a F-law 
EOS with r = 4/3. We perform simulations with uniform 
resolutions Ax = Ay = A = 0.004, 0.0025 and 0.002, ap- 
plying "copy" boundary conditions at the computational 
domain boundaries. We find that evolutions with HLL 
fiux, coupled with either the MC or PPM reconstruction 
schemes, result in a code crash due to the strong initial 
pressure jump. We are able to evolve the system stably 
by using the minmod reconstruction scheme cotipled with 
the LLP fiux. A similar finding is reported in [52|. 

Figure M shows the two dimensional profile of density 
Po, gas pressure P, magnetic pressure b^/2 and magnetic 
field lines at time t — 0.4, where the blast wave has nearly 
reached the boundary of the computational domain. Fig- 
ure [7] shows the one-dimensional profiles along the x and 
y axis for the three resolutions. The profiles are qualita- 
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FIG. 3: Same as Fig.[T]but for the ID shock tube 1 and shock 
tube 2 tests. 



FIG. 5: ID nonhnear Alfven wave test: L2 norms of the errors 
in u^, u^, B^ and B^ at i = fflnai = 2.0. This log- log plot 
demonstrates that L2 norms of the errors are proportional to 
Ax^, and are thus second-order convergent. 




FIG. 4: ID nonlinear Alfven wave test: MHD variable pro- 
files. Test results with resolution Ax = 0.01 (dashed, black 
lines) are compared to the exact solution (solid, red lines) 
at time t — tfinai = 2.0. Our computational domain is 
xG(-2,2). 



tively similar to the those reported in [2 71 . ISll |52| . The 
initial high pressure in the central region causes a strong 
explosion. The explosion is asymmetric in the x— and 
y— directions due to the presence of magnetic fields. The 
explosion is unimpeded in the x-direction, so the Lorentz 
factor u° of the fluid is larger along the x-axis than along 
the y-axis, as demonstrated in Fig. [71 The magnetic field 
lines are squeezed in the y direction, sapping the mag- 
netic field energy in the central region, and driving an 
intense magnetic field in two thin oblate layers surround- 
ing the central region (see Fig. [5]). By t — 0.4, the cen- 
tral density and magnetic pressure have decreased by two 
orders of magnitude, while the central gas pressure has 
dropped by three orders of magnitude. 

By comparing the numerical data from the three reso- 
lution runs in the entire computational domain att = 0.4, 
we see signs of convergence. However, the convergence 
rate is less than first order. This is likely due to the 
fact that the initial strong pressure discontinuity requires 
resolutions higher than those used in our simulations to 
exhibit the proper convergence, as pointed out in [5^ . 
However, we find that in the central \x\ < 0.3, |y| < 0.3 
region, po, P and b^ converge to second order, while u^ 
converges to first order. 

Cylindrical rotating disk (rotor) 

The initial configuration of this rotor test consists of a 
uniform high density (po — 10) central region of cylindri- 
cal radius 0.1 uniformly rotating with an angular velocity 
id = 9.95. The disk is surrounded by an ambient gas of 
density po = 1. The gas pressure P = 1 is constant ev- 
erywhere. The initial magnetic field is uniform and is set 
to S^/V4^ = 1 and Bv = B^ = 0. The gas satisfies 
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FIG. 6: Cylindrical blast explosion: 2D MHD variable profiles. Density po, gas pressure P, magnetic pressure b /2, and 
magnetic field lines are plotted at i = 0.4. The simulation is performed with uniform resolution Aa; — Ay = A = 0.002. 
Magnetic field lines coincide with contours of Az, and are thus plotted according to A^ — 0.5i — 8, with i = 1, 2, . . . , 31. 



a r = 5/3 EOS. We evolve the system using the min- 
mod reconstruction scheme coupled with the LLF flux 
and at resolutions Ax = Ay = A = 0.004, 0.0025 and 
0.002. "Copy" boundary conditions are applied at the 
outer boundaries for this test. 

Figures [5] and [H] show the profiles of po, P, b^f^, u^ and 
magnetic field lines at time t = 0.4. These profiles are 
qualitatively similar to those in [23, [52] • The rotor causes 
magnetic winding. At time t = 0.4, the field lines in the 
central region are rotated by ~ 90°. The winding slows 
down the rotation of the disk. The maximum Lorentz 
factor decreases from the initial value of 10 to 1.7 at 
t = 0.4. The density, pressure and magnetic field in the 
central region also decrease substantially. A high-density, 
oblate shell is formed surrounding the central region. 

As in the cylindrical explosion test, we see signs of 
convergence as the resolution is increased. However, the 
overall convergence rate is less than first order due to 
resolutions too low to adequately resolve the fine struc- 
ture of the flow. The rotor test is even more severe than 



the cylindrical explosion test. This is because the initial 
Lorentz factor u° has a steep slope near the edge of the 
disk. Even with our highest resolution A = 0.002, the 
initial u° decreases from 10 at the edge of the disk to 4.5 
at the next grid point inside the disk. While the three 
simulations produce the same qualitative result, proper 
convergence order is not likely to be achieved when this 
initial steep feature of the velocity is poorly resolved. 
However, we do flnd approximate second-order conver- 
gence in Po and b^ in the the region along the x-axis with 
l^l < 0.2 before the density, pressure and magnetic field 
display a sudden jump (see Fig. [9]). On the other hand, 
M° and P converge faster than first order but less than 
second order in that region. 



There are several conserved global quantities in two- 
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FIG. 7: Cylindrical blast explosion test: ID MHD variable profiles at different resolutions. Density po, pressure P, magnetic 
pressure b^ /2, and Lorentz factor it" profiles along the a;-axis (left) and y-axis (right) at i = 0.4 are plotted at resolutions 
A =0.004 (black solid line), 0.0025 (red dotted line) and 0.002 (blue dashed line). 



dimensional Minkowski spacetime: 



M 



E 



Pk 



J 



Pqu dxdy 
T°"dxdy = M 



E 



p^ijAxAy 






fijAxAy 



T kdxdy 



i,j 



{Sk)ijAxAy , 



{xT°y - yT°x)dxdy , 



(75) 
(76) 
(77) 
(78) 



where the sum is over all the grid points and the volume 
average is equivalent to the surface average over a grid 
cell in the x-y plane in two dimensions. Since there is 
no source term in Minkowski spacetime [i.e. S = in 
Eq. p3l) ]. our finite- volume scheme should conserve Af, 
E, and Pk to roundoff error, provided that no material 



fiows through the boundary of the computation domain 
(i.e. F ^ at the outer boundary). This condition is sat- 
isfied in our rotor test, since the ambient medium is static 
and the torsional Alfven wave generated by the rotor and 
the expansion of the high density gas have not reached 
the boundary at the end of our simulations at t — 0.4. 
Our numerical data confirm that M, E and Pk are indeed 
conserved to roundoff error. On the other hand, the an- 
gular momentum will not be conserved to roundoff error 
since we use Cartesian coordinates to evolve the system. 
Strict numerical conservation of angular momentum can 
be achieved if cylindrical coordinates are adopted (how- 
ever, Px and Py will not be strictly conserved in cylin- 
drical coordinates). We find that for the rotor test at 
t = 0.4, J is changed by 1.7% from its initial value when 
evolved with resolution A = 0.004, 1.2% with A = 0.0025 
and 1.0% with A = 0.002. The slow decrease in J vio- 
lation with resolution is again related to the insufficient 
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FIG. 8: Cylindrical rotating disk (rotor) test: 2D MHD variable profiles. Density po, gas pressure P, magnetic pressure 
6^/2, and magnetic field lines are plotted on the a;y-plane at i = 0.4. The simulation is performed with a uniform resolution 
Aa; — Ay = A = 0.002. Magnetic field lines coincide with contours of Az, and are thus plotted according to Az = 0.16i — 2, 
with i = 1,2, ...,24. 



resolution to resolve the initial steep m° profile near the 
edge of the rotor. We find that the numerically computed 
initial J deviates from the analytic value by 6.8%, 2.7% 
and 1.8% for A =0.004, 0.0025, and 0.002, respectively. 
This indicates that the thin layer near the edge of the ro- 
tor with high initial u^ has a non-negligible contribution 
to J. Angular momentum conservation can be improved 
substantially if the thin layer is well-resolved. 

It follows from the induction equation dtB+'VxE — 
that the global quantities 



B'^dxdy 



(79) 



are conserved as long as E vanishes at the boundary. 
Since we do not evolve the volume-averaged B^, but in- 
stead Eqs. (pD|) - ([S5)) . our scheme does not conserve Q'' 



to roundoff error. Instead, the quantities 



(80) 



are strictly conserved in our scheme. Our numerical data 
confirm this expectation. The deviation between Q'' and 
Q^ converges to zero at second order with increasing res- 
olution. Unlike the angular momentum, the strict con- 
servation of Q'i means that Q^ — Q^ is time independent 
and therefore will not grow with time during the evolu- 
tion. 



B. Curved spacetime test: Relativistic Bondi flow 

Next, we test the ability of our code to accurately 
evolve the relativistic MHD equations in a strongly 
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FIG. 9: Cylindrical rotating disk (rotor) test: ID MHD variable profiles at different resolutions. Density po, pressure P, 
magnetic pressure 6^/2, and Lorentz factor u^ along the x-axis (left) and j/-axis (right) at f = 0.4 are plotted, at resolutions 
A =0.004 (black solid line), 0.0025 (red dotted line) and 0.002 (blue dashed line). 



curved spacetime near a black hole. We perforin the 
magnetized relativistic Bondi accretion test. Bondi ac- 
cretion refers to spherically symmetric, steady-state ac- 
cretion of a unmagnetized, adiabatic gas onto a station- 
ary star. The gas is assumed to be homogeneous and at 
rest far from the star and flow adiabatically with a F- 
law EOS. Analytic solutions for Bondi accretion onto a 
Schwarzschild black hole are given in [5j, [5J| . It has been 
shown that the relativistic Bondi solution is unchanged in 
the presence of a divergenceless radial magnetic field [5a| . 

This test is a powerful one, since it combines strongly 
curved spacetime and relativistic fiows, with an analytic 
solution against which we compare our numerical re- 
sults. It can also be used to test the ability of our AMR 
GRMHD scheme to handle the black hole interior, espe- 
cially the coordinate singularity at the center. The use 
of refinement boxes is natural, since higher resolution is 
required in the vicinity of the black hole, whereas a rela- 



tively low resolution is sufficient to resolve the region far 
away from the black hole. In addition to simulations on 
a fixed background spacetime, we also evolve the black 
hole spacetime using the puncture technique. The spa- 
tial metric then evolves from the puncture initial data to 
the trumpet solution [56, 57]. Although this evolution is 
a pure gauge effect, the spatial metric and extrinsic cur- 
vature change with time, and the gas and magnetic field 
will respond to this change. 

In general, a spherically symmetric spatial metric can 
be written in the form 

'■^Us^ = A{r,t)dr^ + X{r,t)r'^{d9'^ +sm'^ 9d(l)^) . (81) 

It is easy to show that any divergenceless, radial magnetic 
field is given by 



B^ir,t) = 



y/M^ X{r,t)r^ 



(82) 
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where M is the mass of the black hole and Bq is a con- 
stant characterizing the strength of the magnetic field. 
Cartesian coordinates can be constructed from the usual 
transformation: x = r sin 9 cos cj), y = r sin 9 sin (f> and 
z = r cos 9. The Cartesian components of the magnetic 
field B'' is given by 



B'{x,t) 



BaM'^x' 



(83) 



where the determinant 7 of the spatial metric jij in 
Cartesian coordinates is given by ^{x, t) = A(r, i)A^(r, t). 
It is easy to show that this magnetic field can be derived 
from the vector potential 



A^ — 



BgM'y 

r{r + z) 



A, 



r{r + z) 



A,=0 



(84) 



We note that Eq. (|83)) is quite general. The radial coor- 
dinate r can be the Kerr-Schild radius, the shifted Kerr- 
Schild radius considered below, the isotropic radial co- 
ordinate in the puncture initial data, or the radial coor- 
dinate in the trumpet solution of a Schwarzschild black 
hole. During the puncture evolution of a Schwarzschild 
black hole, the radial coordinate r — y^x^ + y'^ + z'^ in 
the numerical simulation changes from the isotropic ra- 
dial coordinate to the radial coordinate of the trumpet 
solution. Equation (j83|) will remain valid if the evolution 
preserves spherical symmetry. 

The parameters of the magnetized Bondi test pre- 
sented here are the same as those used by [5^, [S^, [53| . 
The sonic radius exists at Schwarzschild (areal) radius 
rs = 8M. The density is normalized so that the mass 
accretion rate is Af = 1, and the equation of state is 
r = 4/3. The initial data for the hydrodynamic vari- 
ables are given by the analytic solution, and the magnetic 
vector potential is given by Eq. (1841) . We parametrize 
the strength of the magnetic field by the ratio b'^/po 
at the event horizon. The relationship between Bq and 
(&^/po)horizon cau bc Computed analytically and is given 
by 



Bn = 




r^O / horizon 



(85) 



for the hydrodynamic setup chosen here. We note that 
even though Eq. ((85|) is computed in Kerr-Schild radial 
coordinates, it applies to any other radial coordinate be- 
cause Bq is gauge-invariant. To see this, we compute b^ at 
spatial infinity, where the gas is static. Using Eqs. (|8T|), 
and ITUl) we obtain 



b'^[r) = Jrr 



{B'-f BlM^ 



for r — >■ 00, 



47r 47rr^ 

where ra = vA r is the areal radius and we have used the 
fact that au° = \/l + ■y^^UiUj = 1 for a static {ui = 0) 
gas. Hence we can write 



Bo^ lim 4^ -^ b^ra) 



(86) 



which is manifestly gauge-invariant. 

In all of our simulations, we use five refinement boxes 
with half-side lengths of 3.125M, 6.25M, 12.5M, 25M, 
and 50M. The outermost, lowest-resolution box pos- 
sesses half-side length lOOAf . We only evolve the space 
above and on the equatorial plane z > 0. Equato- 
rial symmetry is applied to hydrodynamic variables and 
(r + z)Ai. All variables at the outer boundary are 
frozen to their initial values. Our standard resolution is 
Ax — Ay = Az = A — 2.5M in the coarsest level. The 
grid spacing A decreases by a factor of two at each suc- 
cessive refinement level, so the resolution on the finest 
level is Amin = M/12.8. For the purposes of testing 
convergence, we also perform a simulation in which the 
resolution is scaled up so that Amin — M/16. To mea- 
sure errors due to moving refinement boxes in our AMR 
scheme, we move refinement box centers according to 



Xc = a;,„ sinwi , j/c = 2/m(l - coswi) 



(87) 



where we set the parameters Xm = l.OA/, y„i = 0.6Af 
and 2tt/uj — 50M. Below, we present results for the 
fixed background spacetime simulations and the puncture 
evolution. Without loss of generality, we set M = 1 in 
all of our simulations. 



1. Fixed background spacetime 

In many relativistic Bondi tests, Kerr-Schild coordi- 
nates are used together with excision. A different ap- 
proach is adopted here. We first define a shifted Kerr- 
Schild radius r = vks — ^0, where vks is the Kerr- 
Schild radius and tq is a constant chosen in the range 
< ro < 2M. We then construct Cartesian coordinates 
using the standard transformation between (x, y, z) and 
(r, 9, (/)). The origin x = y = z = therefore corresponds 
to a Kerr-Schild radius rxs = '''o- The region tks < ro 
is excluded in this coordinate system and so is the black 
hole spacetime singularity. However, the origin is a co- 
ordinate singularity since the whole surface rKs = '"0 is 
mapped to a single point. This coordinate system thus 
mimics the trumpet solution of a Schwarzschild black 
hole. We set ro = M for all the tests presented in this 
section. Just as in puncture evolutions, we do not use 
excision but shift the grid slightly so that the origin is 
not on a grid point. 

We are able to evolve the system stably using the MC 
reconstruction scheme coupled with the HLL flux for 
(6^/po)horizon S 10 with our Standard resolution. Higher 
magnetic fields may be evolved if the resolution is in- 
creased. During the evolution, the magnetic field, as well 
as the density, increases linearly with time near the ori- 
gin, creating jumps in the magnetic field that increase 
with time. This phenomenon causes the evolution near 
the origin to become more and more inaccurate. The in- 
accurate data spread out slowly from the origin to the 
apparent horizon, eventually crossing into the BH exte- 
rior. To overcome this difficulty, we add fourth-order 
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FIG. 10: Fixed-background, (&^/po)hor: 



=4 magnetized 



Bondi test: ID MHD variable profiles, po, v^ and B^ are plot- 
ted in the equatorial plane (2 — 0) along the line y = O.OlAf 
at t = 101. 25M. Solid (red) lines are the analytic solution and 
dashed (black) lines are numerical data with Amin = M/12.8 
on the finest refinement level, using MC reconstruction. The 
vertical lines denote the location of the event horizon I a:: I — M. 
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FIG. 12: Fixed-background magnetized Bondi test: Conver- 
gence study using MC reconstruction. L2 norm of B^ as a 
function of (6^/po)horizoii is plotted at i = 101.257\/ for lower 
(Amin = M/12.8) and higher (Amin = M/16) resolution runs. 
The lower resolution result is multiplied by the factor 0.64 to 
demonstrate second-order convergence. 
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FIG. 11: Same as Fig. [TOl but with PPM reconstruction for 
the numerical data, and setting Ai = Q for deep inside the 
BH (r < 0.5M). 



Krciss-Oliger dissipation to the evolved variables inside 
the horizon for radius r < Q.8M. We also set a density 
cap Po < 1 for radius r < 0.5Af . This technique stabilizes 
the evolution near the origin and the system quickly set- 
tles down to a steady state inside the horizon. Figure [TU] 
shows the profiles of po, w^ and B^ in the equatorial plane 
{z = 0) along the line y — 0.01 for (6^/po)horizon =4 at 
t = 101. 25M, by which time the center of the refinement 
boxes has gone through slightly more than two rotations. 
The analytic solution and numerical data are plotted to- 
gether for comparison. Vertical lines denote the location 
of the horizon |x| = M. We see that the profiles agree 
very well with the analytic solution outside the horizon. 
There are strong jumps in the magnitude and direction 
of the magnetic field near the coordinate singularity at 
the origin. The maximum and minimum values of B^ 
near the origin are 13.4 and -19.5 respectively, far out- 
side the scale shown in the figure. However, these jumps 
are always contained near the coordinate singularity. 

The evolution near the coordinate singularity at the 
origin is less stable when evolved with PPM reconstruc- 
tion. To remedy this, we set ^^ = for radius r < 0.5M, 
well inside the horizon, in addition to the technique de- 
scribed above. Figure [TT] shows profiles of MHD variables 
using PPM. We again see that the profiles agree well with 
the analytic solution outside the horizon, oblivious to the 
ruggedness of profiles in the black hole interior. 
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To check for convergence, we perform a number of sim- 
ulatfons of varying (6^//Co)horizon with two different reso- 
lutions Amin = M/12.8 and Amin ~ A//16. We compute 
the L2 norm of B"^ a.t t = 101. 25M by summing over grid 
points 



L2(S'') = 



V numerical analytic 



EB, 



analytic,! 



where Bl 



-"analytic,! dcuotc the analytic values of B^ for 
(&^/po)horizon =1- Wc Only computcd the L2 norm in 
the innermost refinement level outside the horizon with 
|a;| < 3M, \y\ < 3M, 0< z < 3M and r > M. This is the 
region in the black hole exterior where the magnetic field 
is the strongest. Figure [T^ shows the L2 norm as a func- 
tion of (&^/po)iiorizon • Wc sce sccoud-ordcr convergence 
for (6^//3o)horizon ^ 8. The convergence rate appears to 
be higher than second-order for (6^/po)horizon =10, indi- 
cating that the data in the lower resolution run may not 
be accurate enough to display proper convergence. 



2. Puncture evolution 
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In addition to the evolution with a fixed background 
metric, we perforin several magnetized Bondi tests with a 
time- dependent background metric. We evolve the black 
hole spacetime using the puncture technique. In order 
to compare with the analytic solution, we set the matter 
and EM field source terms to zero in the BSSN equations, 
so that the gas and EM field do not affect the spacetime 
evolution, consistent with assumptions used when deriv- 
ing the analytic solution. However, the gas and EM field 
will respond to the change of the background metric since 
the MHD and induction equations contain metric quan- 
tities. The metric evolution is a pure gauge effect: the 
spatial metric and extrinsic curvature evolve from the 
initial maximal, wormhole slicing to the final slicing rep- 
resenting the trumpet geometry. 

We evolve the MHD and induction equations using 
both MC and PPM reconstruction, coupled with the HLL 
flux. A fourth-order Kreiss-Oliger dissipation is applied 
to the MHD evolution variables for r < 0.5 M, which is 
inside the horizon at all times. As before, we set Ai = 
for r < 0.5M in the PPM run to stabilize the evolution 
near the puncture. Figure [T3] shows the profile of B'^ at 
t — 101. 25M in the equatorial plane {z = 0) along the 
line y = O.OIM for (&^/po)horizon =4. Numerical data 
are compared to Eq. ((55)) with y/j — e^"^ taken from the 
numerical data. We see that the data agree well with 
the analytic result outside the horizon in both runs. The 
glitch near the origin results from the loss of accuracy of 
(j) near the puncture. When compared with Fig. IIOI we 
see that the B^ profile is smooth in the puncture evolu- 
tion with MC reconstruction. However, we find a similar 
feature in the po profile as in Fig. [TOl 

Since the analytic solution of the hydrodynamic quan- 
tities are given in Kerr-Schild coordinates, direct compar- 
ison of numerical and analytic results is not easy in these 



FIG. 13: Evolved-spacetime, (6 /po)horizon =4 magnetized 
Bondi test: Profile of B"" at i = 101. 25M in the equato- 
rial plane {z = 0) along the line y = O.OIM. The metric 
is evolved using the puncture technique. The upper graph 
plots numerical data using MC reconstruction, and the lower 
graph shows the result using PPM reconstruction and setting 
Ai — for r < 0.5M. Dashed (black) lines are numerical 
data, and solid (red) lines are results computed by Eq. (|83p 
with y/^ = e^'^ taken from numerical data. The glitch near 
X = results from the loss of accuracy of the metric data 
close to the puncture. The vertical lines denote the location 
of the black hole horizon la;! = 0.94M. 



simulations. However, since both b^ and po a-re scalar and 
the system is stationary, the profile of 6^ as a function of 
Pq is gauge-independent. Figure [HI shows this function 
at i = 101. 25M in the equatorial plane along the line 
y — 0.01 and x > 0.5M for two resolutions. The numeri- 
cal profile of po is no longer monotonically increasing with 
decreasing r when the numerical data inside the horizon 
are included, due to inaccuracy near the puncture. We 
therefore remove the data points inside the horizon for 
x < xo to prevent multiple values of 6^(po) from appear- 
ing in the plot, where xq is the point when po reaches the 
maximum. The position of the horizon is indicated by 
the vertical line po = (po)horizon = 0.02579, the value of 
po at the horizon. The deviation between the numerical 
data and analytic result becomes visible close to the hori- 
zon in the lower resolution run Amin = M/12.8. Much 
better agreement is achieved in the higher resolution run 
with Amin = M/16. This is not surprising since Af/12.8 
is a fairly poor resolution for puncture simulations. 
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FIG. 14: Evolved-spacetime, (6 /po)horizon =4 magnetized 
Bondi test: Convergence of \? as a function of po at t = 
101.25Af along the line x > 0.5M, (j/, z) = (0,0.01M). The 
background metric is evolved using the puncture technique. 
Dotted (black) and Dashed (blue) lines plot the numeri- 
cal data evolved with MC reconstruction using resolutions 
Amin = M/12.8 and Amin = M/16, respectively. The solid 
(red) line denotes the analytic profile, and the vertical line 
demarcates the horizon boundary, where po = (po) horizon ~ 
0.02579. The region with po > (po)horizon lies inside the hori- 
zon. 



C. Curved spacetime test: Collapse of magnetized 
rotating relativistic star 

This test focuses on magnetized, rotating, relativistic 
stellar-collapse simulations. The initial stellar configura- 
tion is the same as Star D in W] and Star B in |1'] . The 
star satisfies a F = 2 polytropic EOS and is uniformly 
rotating with J/AP = 0.34, where J is the angular mo- 
mentum. The ADM mass of the star is M = 1.04Mtov, 
where Mrov is the maximum ADM mass of a non- 
rotating relativistic star satisfying F = 2 EOS. The star 
is on the unstable branch of the constant J sequence, and 
previous numerical simulations have demonstrated that 
it is dynamically unstable to gravitational collapse [l|,l60|. 

In all of our simulations, we use seven refinement boxes 
with half-side lengths of 0.9143M, 1.829M, 3.657M, 
7.314M, 14.63M, 29.26M and 58.51Af . The box contain- 
ing the outer boundary has half-length 117.0Af . The ini- 
tial coordinate radius of the star in the equatorial plane is 
3.485M. Hence the stellar interior is initially covered by 
the three innermost refinement boxes. The grid spacing 
is reduced by a factor of two at each successive refinement 
level. We perform three simulations with the resolution 



FIG. 15: Magnetized stellar collapse test: Evolution of the 
central lapse Oc for the low (black solid line), medium (red 
dotted line), and high (blue dashed line) resolution runs. The 
dot in each case indicates the time at which the apparent 
horizon appears. The increase in ac soon after the horizon 
formation is caused by the loss of accuracy in metric evolution 
near the newly formed puncture, which is located near the 
coordinate origin and is deep inside the horizon. 



in the finest refinement level set to Amin = 0.02857M 
(low resolution run), 0.02287Af (medium resolution run), 
and 0.01829Af (high resolution run). 

We evolve the metric using a fourth-order finite- 
differencing scheme. We adopt the puncture gauge condi- 
tions with the shift parameter rj set to 0.5/Af . The MHD 
and induction equations are evolved using the PPM re- 
construction scheme coupled with the HLL flux. Equa- 
torial symmetry is applied to all variables. We maintain 
a low density atmosphere in the computational domain 
with patm = 10 ~^°Pnia x(0) and Patm = -Pcoid(patm) as de- 
scribed in Sec. IIIIEl The Sommerfeld outgoing wave 
boundary condition is applied to the BSSN evolution 
variables, and outflow boundary conditions are applied 
to the hydrodynamic primitive variables, while the vec- 
tor potential Ai is linearly extrapolated to the boundary. 

Since the star is unstable, collapse can be triggered by 
numerical truncation error during the evolution. How- 
ever, since the truncation error is reduced with increas- 
ing resolution, subsequent evolution of the star will de- 
pend sensitively on resolution, which is not desirable for 
a convergence test. We therefore induce the collapse by 
depleting the initial pressure by one percent. We set up 
a small, poloidal, axisymmetric magnetic field by setting 
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FIG. 16: Magnetized stellar collapse test: Evolution of the 
irreducible mass Mi^, black hole mass Afan and angular mo- 
mentum Jbh, as normalized by the initial ADM mass Af and 
angular momentum J. Shown here are data from the high 
resolution run. Results from the low and medium resolution 
runs are similar. 



the vector potential as follows: 



Ax = -y Ab ma,x{P ~ Pcut,0) 
Ay = X Ab vaeoi{P - Pcut.O) , 
A, ^ 0, 



(89) 
(90) 
(91) 



where Pcut is set to 4% of the initial maximum pressure. 
The constant parameter Ab determines the strength of 
the magnetic field. We characterize the strength of the 
magnetic field by the ratio of the magnetic energy A4 to 
the internal energy i?int- These energies are defined as 



Eint = / y/^{poe)u d X , 



(92) 
(93) 



We have chosen a magnetic field strength of M/Ei^t = 
7.3 X 10~^, which introduces only a small perturbation 
to the star. 

Figure [15] shows the evolution of the central lapse for 
the three resolution runs. As the star is collapsing, the 
lapse decreases and an apparent horizon appears at time 
t ^ TOM. The large energy density and magnetic pres- 
sure inside the horizon causes the code to crash soon after 
its formation. This difficulty can be overcome by evacu- 
ating the hydrodynamic matter and magnetic field deep 
inside the horizon soon after the formation of horizon. 



FIG. 17: Magnetized stellar collapse convergence tests. Up- 
per graph: Magnetic pressure Pmag = fe^/2 as a function of 
X along the diagonal line x = y = z at time t = 40.2M 
for the low (black solid line), medium (red dotted line) and 
high (blue dashed line) resolution runs, normalized by the 
initial maximum value of Pmag- Lower graph: Pairwise differ- 



ences of Pmag between different resolution runs. The differ- 
ence SP^'-'^^'' = (P^k^ - P,S'at)/^mag,max(0) is multiplied by 
1.5625 to demonstrate deviations from second-order conver- 
gence. Notice that the results converge slightly higher than 
second order in the high Pmag region but less than second 
order in the low Pmag region. 



The evolution then proceeds stably, and the spacetime 
settles to a Kerr black hole after t > 75M (see Fig. [T6|) . 
with virtually no fiuid or magnetic fields left outside the 
horizon. The mass and spin of the black hole are com- 
puted using the isolated and dynamical horizon formal- 
ism [6l| , with the axial Killing vector field computed us- 
ing the numerical technique described in [62]. We find 
A/bh ~ M and Jbh « J (asH/MBH = J/M"^ = 0.34) 
for all three resolution runs once all the matter enters 
the horizon, where M and J are the initial ADM mass 
and angular momentum of the star, respectively. This is 
expected since the collapse is nearly axisymmetric, and 
only a negligible amount of mass as well as angular mo- 
mentum is radiated by the gravitational waves. (Recall 
that no angular momentum is radiated in strict axisym- 
metry.) 

Figure [17] shows the profile of magnetic pressure 
^mag = b^/2 along the diagonal line x = y = z ai 
t = A0.2M. We see slightly higher than second-order con- 
vergence in the high Pmag region but lower than second- 
order in the low Pmag region. 

In the simulation with the highest resolution, we find 
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that the vector potential Ai develop spikes near the sec- 
ond innermost refinement boundary during and after the 
edge of the Ai — surface passes through that refine- 
ment boundary. The amplitude of the spikes amplifies 
with time, eventually causing the code to crash. This dif- 
ficulty can be removed by adding a fourth-order Kreiss- 
Oliger dissipation to Ai . The origin of the spikes is from 
prolongation and restriction. As Ai are steeply decreas- 
ing to zero near the edge, our adopted third-order La- 
grangian interpolation scheme adds spurious oscillations 
in Ai near the refinement boundary after prolongation 
and restriction. Since the refinement boxes are not mov- 
ing, the oscillation amplitude amplifies each time when 
prolongation and restriction are applied. The same phe- 
nomenon could occur for other hydrodynamical variables 
with a steep gradient. However, this effect has a more sig- 
nificant impact on the magnetic field, since a slight spa- 
tial oscillation in Ai will be amplified after taking spatial 
derivatives. An alternative method to cure this prob- 
lem would be to use a more sophisticated interpolation 
scheme such as the ENO or WENO scheme. We plan to 
investigate these alternative interpolation schemes in the 
future. 



V. CONCLUSION 

We have developed a new GRMHD code that is capable 
of evolving MHD fiuids in dynamical spacetimes. We 
use the BSSN scheme coupled with the puncture gauge 
conditions to evolve the metric, and an HRSC scheme to 
evolve the MHD and induction equations. 

We adopt the formalism described in [27] to recast the 
induction equation into an evolution equation for the 
magnetic vector potential Ai [i.e. Eq. (IT^ ]. The vari- 
ables Ai are stored on a staggered grid with respect to the 
other variables. The divergenceless constraint 'V ■ B — 
is imposed through the vector potential. This evolution 
scheme is AMR-compatible, with prolongation and re- 
striction applied to the unconstrained variables Ai in- 
stead of B*, which gives us fiexibility in choosing dif- 
ferent interpolation schemes for prolongation/restriction. 
In simulations with uniform grid spacing, our scheme for 
evolving the magnetic field is numerically equivalent to 
the commonly used constrained-transport scheme based 
on a staggered mesh algorithm |16| . 

We have performed several code tests to validate our 
code, including magnetized shocks, nonlinear Alfven 
waves, cylindrical blast explosions, cylindrical rotating 
disks, magnetized Bondi tests, and collapse of magne- 
tized rotating stars. We find good agreement between 
the analytic and numerical solutions, and achieve second- 
order convergence for smooth flows, as expected. 

In GRMHD simulations in dynamical spacetimes in- 
volving black holes, one delicate issue is the handling of 



the black hole interior. We adopt the moving puncture 
technique in which the black hole spacetime singularity 
is avoided by the puncture gauge conditions. However, 
a coordinate singularity (puncture) remains in the black 
hole interior, which could cause numerical difficulties in 
MHD simulations. In our tests involving black holes, we 
find that the evolution in the black hole interior is more 
stable when a more diffusive scheme such as the MC re- 
construction scheme is used rather than the PPM scheme. 
We plan to investigate the idea of using a less diffusive 
scheme (such as PPM reconstruction coupled with the 
HLL flux) in the black hole exterior and a more diffusive 
scheme (such as MC or minmod reconstruction coupled 
with the LLF flux) in the black hole interior. A similar 
technique is used in some MHD simulations of magne- 
tized accretion disks around a black hole J)3j. We also 
find that adding Kreiss-Oliger dissipation to MHD vari- 
ables in the black hole interior can stabilize the evolution. 

In GRMHD simulations using an FMR grid, we find 
that applying a high order interpolation scheme on Ai 
during prolongation and restriction could cause oscilla- 
tions in Ai near the refinement boundaries. The oscilla- 
tion amplitude can amplify with time. This numerical ar- 
tifact degrades the accuracy of the simulation and could 
even cause the code to crash. The artifact can be re- 
moved by adding a fourth order Kreiss-Oliger dissipation 
to Ai. A better solution is to use a more sophisticated 
interpolation scheme for Ai , such as the ENO or WENO 
scheme. We plan to investigate these alternative schemes 
in the future. 

In addition to the treatment of the black hole interior, 
our MHD code has limitations similar to those of other 
MHD codes in the literature. In particular, accurate evo- 
lution is difficult when b^ 3> po- This could potentially 
cause problems in the low-density regions in some appli- 
cations. However, our experience and the experience of 
other numerical MHD groups suggests that these difficul- 
ties are surmountable. 

Having demonstrated the validity of our AMR 
GRMHD code, we will next apply our code to study 
the effects of magnetic fields in the coalescence of bi- 
nary neutron star and black hole-neutron star systems, 
the collapse of magnetized supermassive stars, and the 
dynamics of magnetized accretion disks around merging 
binary black holes. 
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